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ABSTRACT 

We consider a self-similar force-free wind flowing out of an infinitely thin disk 
located in the equatorial plane. On the disk plane, we assume that the magnetic stream 
function P scales as P oc B? , where R is the cylindrical radius. We also assume that the 
azimuthal velocity in the disk is constant: v$ — Mc, where M < 1 is a constant. For 
each choice of the parameters v and M, we find an infinite number of solutions that are 
physically well-behaved and have fluid velocity ^ c throughout the domain of interest. 
Among these solutions, we show via physical arguments and time-dependent numerical 
simulations that the minimum-torque solution, i.e., the solution with the smallest 
amount of toroidal field, is the one picked by a real system. For v ^ 1, the Lorentz 
factor of the outflow increases along a field line as 7 ~ M(z/R{ p )^ 2 ~ 1 '^ 2 « R/Ra, 
where Rf p is the radius of the foot-point of the field line on the disk and Ra — Rf P /M 
is the cylindrical radius at which the field line crosses the Alfven surface or the light 
cylinder. For v < 1, the Lorentz factor follows the same scaling for z/Rf p < M^ 1 '^ 1- "' , 
but at larger distances it grows more slowly: 7 w (z / Rf p ) u / 2 . For either regime of i>, 
the dependence of 7 on M shows that the rotation of the disk plays a strong role in jet 
acceleration. On the other hand, the poloidal shape of a field line is given by z/Rf p « 
(R/R{ p ) 2 /( 2 ~^ and is independent of M. Thus rotation has neither a collimating nor 
a decollimating effect on field lines, suggesting that relativistic astrophysical jets are 
not collimated by the rotational winding up of the magnetic field. 
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1 INTRODUCTION 

Although relativistic jets from accreting black holes have 
been known and have been studied for decades, the physics 
behind the production, acceleration and collimation of these 
jets is still poorly understood. Broadly, there are two schools 
of thought on the launching of jets. According to the first 
(e.g., Lovelace 1976), jets are the innermost and most en- 
ergetic region of an extended outflow from the black hole 
accretion disk. The power for the jet thus comes from the 
disk. In the other view, the jet is powered by the free energy 
associated with the spin of the central black hole (Blandford 
& Znajek 1977). The jet phenomenon is then a remarkable 
manifestation of the Penrose (1969) process. 

The fact that jets from black holes are observed to ac- 
celerate to relativistic velocities, despite the absence of any 
significant radiative or thermal driving, strongly suggests 
that magnetic fields play a dynamically important role. This 
motivates the study of magnetohydrodynamic (MHD) mod- 
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els of jets. Moreover, given the difficulty of achieving self- 
collimation in an isolated jet — e.g., Michel's (1973a) solu- 
tion for a spinning monopole has no collimation at all (see 
§ 3.2 below) — it is reasonable to suppose that the pres- 
ence of an extended outflow surrounding the jet is needed 
for maintaining both the collimation and dynamical stabil- 
ity of the jet, as stressed for instance by Appl & Camenzind 
(1992, 1993) and Beskin & Malyshkin (2000). This external 
flow is presumably a wind of some sort from the disk. It is 
thus of interest to study the dynamics of extended, nearly 
self-similar, MHD, disk outflows. 

Given the complexity of the full relativistic MHD equa- 
tions, a variety of simplifications have been introduced by 
different authors. Many studies make a nonrelativistic ap- 
proximation (e.g., Blandford & Payne 1982; Heyvaerts & 
Norman 1989; Contopoulos & Lovelace 1994; Contopoulos 
1995b; Ostriker 1997). Although it is adequate for certain 
aspects of the problem, this approach fails to capture the 
presence of a light cylinder (or relativistic Alfven surface), 
an explicitly relativistic phenomenon that likely plays an im- 
portant role in jet physics. The reader is referred to the semi- 
nal work on light cylinders by Goldreich & Julian (1969) who 
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discussed pulsar magneto-spheres and Okamoto (1974, 1978) 
who considered disk winds, as well as the work of Lovelace 
et al. (1986), Nitta, Takahashi & Tomimatsu (1991), Beskin 
& Pariev (1993) and Beskin (1997) who included the effects 
of general relativity. 

In considering relativistic winds, one could either work 
with the full MHD equations (e.g., Vlahakis & Konigl 2003), 
or make further simplifications. One possibility is to assume 
that the outflowing gas is cold (e.g., Okamoto 1978; Li et 
al. 1992; Contopoulos 1994) so that only the inertia of the 
fluid contributes to the dynamics and there are no effects 
due to pressure or internal energy. An even more drastic 
simplification is to ignore gas inertia as well. In this force-free 
approximation, the plasma supplies charges and currents as 
needed to support the electromagnetic fields, but it has no 
other dynamical role. 

Force-free electrodynamics is the simplest and clean- 
est model of magnetized relativistic outflows. If magnetic 
fields are at all important in disk winds and jets, one might 
hope that the force-free approximation would capture a good 
fraction of the relevant physics. Early discussions of force- 
free disk winds may be found in Okamoto (1974), Blandford 
(1976) and Blandford & Znajek (1977). Recently, force-free 
electrodynamics has been formulated as a syst e m of time- 
dependent equatio ns jKomissarovll2002bL l2004t iMcKinnevI 
2006a| ISpitkovskvl [20061. These time-dependent force-free 
equations have been use d to study black hole a nd neu- 
tron st ar magnetospheres (|Kom issarovll200lll20 02allbl,l2004 |. 
| 2006l: IMcKinnevI 12006 alibi llVIcKinnev fc Naravanl l2006bl: 
ISpitkovskvll2006h . 

Force-free models have gained new relevance as a result 
of time-dependent general-relativistic MHD simulations of 
black hole accretion flows (e.g., McKinney & Gammie 2004; 
Gammie et al. 2004; de Villiers et al. 2005). These simula- 
tions reveal a nearly force-free highly relativistic jet form- 
ing spontaneously along the symmetry axis and accelerating 
rapidly outwards. Moreover, the height-integrated azimuthal 
current within the turbulent disk is found to have a sur- 
prisingly simple power-law behavior (McKinney & Narayan 
2006a, b), dl^/dR oc R~ 5 ^ 4 , where R is the radius and I^iR) 
is the azimuthal current enclosed inside R. In addition, an 
idealized force-free model with an equatorial rotating cur- 
rent sheet with azimuthal current varying as P _5/4 , which 
is equivalent to the magnetic stream function varying as 
P oc i? 3,/4 (or v = 3/4 in eg. 1141 below ), gives a surprisingly 
good match to a variety of properties of the simulated jet 
such as the poloidal configuration of the magnetic field, the 
collimation angle of the jet, the Lorentz factor, etc. Thus, 
it appears that simple disk wind models with self-similar 
(power- law) scalings may capture key aspects of the jet prob- 
lem. We should note, however, that self-similar disk wind 
models with the particular scaling P oc P 3//4 have a check- 
ered history and there has been some confusion on whether 
or not well-behaved solutions that extend to large distances 
from the disk are present at all. Blandford & Payne (1982) 
and Contopoulos (1995a), for instance, both find that such 
solutions do not exist. 

In this paper we analyze the structure of self-similar 
force-free disk winds for a range of power-law profiles of the 
stream function, including the case P oc R 3/4 . § 2 introduces 
our model, which is closely similar to the model of Contopou- 
los (1995a) and is a self-similar specialization of Okamoto's 



(1974) more general analysis. We show that the problem 
is reduced to solving the second-order differential equation 
126H with appropriate boundary conditions. § 3 discusses a 
number of analytical results and maps out the different kinds 
of solutions that are possible. § 4 presents results from a nu- 
merical analysis of the basic differential equation and iden- 
tifies regions in parameter space where the different solution 
types are found. § 5 presents time-dependent numerical sim- 
ulations of force- free disk winds using a relativistic force- free 
code and draws some general conclusions as to which of the 
many solutions are relevant for given situations. The paper 
concludes with a discussion in § 6. The Appendix discusses 
some mathematical properties of the differential equation 



2 THE MODEL 

2.1 Okamoto's (1974) Force-Free Field Equation 

We consider a steady, axisymmetric force-free field in cylin- 
drical coordinates: R, <f), z. Since V • B = 0, the mag- 
netic field B can be written in terms of a vector potential: 
B = V x A. Furthermore, since we have assumed axisymme- 
try (d/d(p = 0), the poloidal component B p of the magnetic 
field (the projection on the Rz plane) can depend only on 
the <fi component of the vector potential. Thus, we may 
write the field quite generally as 



B = B p + B^cj) = V x (A^) + B^ 

1 d(RA^) 1 d{RA^) 

ti, 



R dz 



0i 



R dR 



(1) 



where B^ is the toroidal component of the field. It is stan- 
dard to rewrite this in terms of the stream function P = 
RAtf, (the stream function is also often represented by the 
symbol if)): 

1 dP „ 1 dP] 



B 



R dz ' 



B 



0i 



RdR 



(2) 



The magnetic flux enclosed within radius R is given by 



$(P) = / B z 2-KR'dR' = 2irP{R). 



(3) 



Thus, apart from a factor of 2ir, the stream function P is 
the same as the enclosed magnetic flux. 

Since the enclosed flux interior to a given field line is 
constant, independent of the z at which it is measured, it is 
clear that P must be a constant along each field line. This 
can also be shown explicitly from equation @ by verifying 
that VP • B = 0. Thus, each field line is uniquely specified 
by its value of P. Another quantity that is constant on each 
field line is the angular velocity fi. Since we imagine that 
the field line is embedded in an accretion disk and corotates 
with it at its foot-point, the angular velocity of the entire 
field line out to infinity is determined by the rotation of the 
disk at the foot-point. Q. is clearly a function only of P. 

By definition, in the frame rotating with the angular 
velocity Q of a field line we are in the comoving frame of the 
fluid and hence there is no electric field (because we assume 
infinite conductivity). Therefore, back in the nonrotating 
frame, there must be an electric field, which is given by 

E = X B = — {-B„ 0, B R ) = ——VP. (4) 

c c c 
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The electric charge density may be obtained from Gauss's 
law: 



n_ 

Aire 



E 



(5) 



d 2 P 



1 dP 



d 2 P 



OR 2 RdR dz 2 



4nc ] 1 dP 



In steady state the current density is given by 
c 



J = 3p+34><i 



47T 



VxB, 



(6) 



where we have once again decomposed the vector into 
poloidal and toroidal components. This gives 

c ^ c fd 2 P 1 dP d 2 P\ , ns 

H = ^xBU = -—^—-- M + —)(7) 

In the force-free limit we have no inertial terms in the 
momentum equation, so the sum of all the electromagnetic 
forces at each point must be zero. Thus the equation of mo- 
tion takes the simple form 



p e E + - j x B = 0. 

c 



(9) 



Since E has no toroidal component, neither can j x B. Thus 
j p and B p have to be parallel to each other. Comparing the 
expression for j p in equation with that for B in equation 
we see that the quantity RB^ must be a unique function 
of P. Hence this is a third quantity that is constant along a 
field line: 



RB d 



c dp 



Jp ~ ^dP Bp - 



(10) 
(11) 



Note that, at each point, j3 is proportional to the net en- 
closed current in the z-direction 

Taking the dot product of equation @ with E, we ob- 
tain 



p e E 2 + -E-(jxB). 



(12) 



Substituting the expressions written down earlier for the 
various quantities, and after some algebra, we finally obtain 
the following partial differential equation for P: 



= 



1 - 



fl 2 R 2 \ ir_p 

or 2 



1 + 



n 2 R 2 \ 1 dP 
RdR 



_ tfR^\ d 2 P dp _ VR 2 dQ 3 
1 ' 1 dz 2p dP c 2 dP ] 1 • 



(13) 



This is the force-free field equation of Okamoto (1974, com- 
pare with his eq. 43; simpler versions of this equation with- 
out the last term were discussed earlier by Mestel 1973, 
Michel 1973b and Scharlemann & Wagoner 1973). 

2.2 Self-Similar Disk Wind 

The discussion so far was quite general. All we assumed were 
the following conditions: (i) steady state, (ii) axisymmetry, 
(iii) force-free (i.e., negligible inertia), (iv) infinite conduc- 
tivity. Now we explicitly consider the problem at hand, viz., 
a self-similar force-free wind flowing out of an accretion disk 
(Contopoulos 1995a). We assume that the disk is infinitely 



thin and located at the equatorial plane, z = 0. Each field 
line is identified by the radius R{ p of its foot-point on the 
disk. The disk supplies two boundary conditions at the foot- 
point: 

1. The magnetic flux &(R{ P ) enclosed inside radius Rf p : This 
determines the stream function P(R{ P ) = Q(Rf p ) /2n of each 
field line as a function of the radius of its footpoint. 

2. The angular velocity Q(Rf p ) of the disk at radius R{ p : The 
foot-point of the field line is dragged around by the disk at 
the angular velocity £l(Rf p ) and since Q is constant along 
each field line this angular velocity is imparted to the entire 
line out to infinity. 

Assuming that we have a self-similar disk with a power- 
law structure, we write the stream function on the disk plane 
as 



P(R,z) 



oc R v , 



s£ v C 2. 



(14) 



By this definition, v is equivalent to the index x in Con- 
topoulos (1995a) and the index F in Vlahakis & Konigl 
(2003). The vertical component of the magnetic field on the 
disk plane then scales as 



B z (R,z) 



1 dP 
R~dR 



oc R"- 



(15) 



Since we would like the enclosed magnetic flux to vanish 
as R — ► 0, we impose the condition v ^ 0. Also, since we 
are not interested in solutions in which the magnetic field 
strength increases with increasing radius, we restrict our- 
selves to v ^ 2. Note that the monopole solution of Michel 
(1973a) corresponds to v = 0, the problem considered by 
Blandford & Payne (1982) corresponds to v — 3/4, the 
paraboloidal field model of Blandford (1976) corresponds 
to v = 1, and the models discussed by Ostriker (1997) cor- 
respond to the range 1 < v < 2. 

The power-law form of the boundary condition I14H mo- 
tivates us to consider the following self-similar form for the 
stream function, 



P(R,z) = R v T{u), 



(16) 



where the self-similar coordinate u (called Z by Contopoulos 
1995a) is given by 



z 

U =R- 

For convenience we set 
T(0) = 1, 



(17) 



(18) 



which is equivalent to absorbing any coefficient on the right- 
hand side of equation 1161 into the definition of R. Since field 
lines live on surfaces of constant P, a field line going through 
the point (7?, z) has its foot-point at radius 



Rf p = RT 1,v {z/R) 



(19) 



Equivalently, for a given Rf p , the shape of the field line in 
the poloidal plane is given by the parametric equations 

R(u) = R Cp T' 1/v (u), z(u) = R ip uT- 1/v {u). (20) 

Consider now the angular velocity. From equation 113H . 
it is clear that the dimensionless ratio QR/c plays a promi- 
nent role in this model. In order to obtain a self-similar 
solution, this quantity has to be constant on the disk plane. 
Thus we write 
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= M = constant, 



(21) 



where M (= l/c n in Contopoulos 1995a) represents a sort 
of "Mach number" of the angular motion. Note that the 
azinmthal velocity follows a "fiat rotation curve," not a Ke- 
plerian profile. This is not the ideal choice for an accretion 
disk, but it is the only profile that is consistent with self- 
similarity (e.g., Li et al. 1992; Contopoulos 1995a; Vlahakis 
& Konigl 2003) . Since Q is constant on a field line, the 
angular velocity at any general (R, z) is given by 



OR 

c 



MT 



-i/ v 



(22) 



Consider now the third conserved quantity on a field 
line, the quantity j3. Under self-similarity, we expect the ra- 
tio of the toroidal and poloidal field strengths at the disk 
plane to be a constant. Thus we assume that 



B,t, 



-sM = constant. 



(23) 



We include M on the right-hand side since we expect the 
amount of toroidal field in the wind to be roughly propor- 
tional to the rotation rate of the disk. We also introduce a 
negative sign because we expect field lines to be swept back 
with respect to the direction of rotation. (The quantity Ho 
in Contopoulos 1995a is —sM in our notation.) We then 
obtain 



P 



P 
dJ3_ 

IF 



-vsMR v ~"T 



v — \rp{y — V) I v 



= v{y-\)s A M z R 



2 jjV — 2rp(v — 2)/v 



(24) 
(25) 



Substituting the various self-similar scalings written 
down above into equation 1131 . we finally obtain the dif- 
ferential equation satisfied by the similarity function T(u): 







(it 2 + 1)T" + (3 - 2v)uT' + v{v - 2)T 



-(v + 2)/v 



(u 2 + 1)TT" - (u 2 + l)T u /v 



+(3 - 2u)uTT' + v(y - 1)(1 - s 2 )T 2 



(26) 



This is the fundamental equation that we need to solve to 
study self-similar force-free disk winds. Once we have a solu- 
tion for T(u), we can calculate all other quantities of interest. 
For instance, the three components of the magnetic field are 
given by 



B R (R,z) 
B4R,z) 
B z (R,z) 



-R v - 2 T'(u), 



-vsMR v ~'T 



v — 2rp{y — \ ) J v 



(«)> 



R v - 2 [vT{u)-uT'(u)], 



(27) 
(28) 
(29) 



where as always it = z/R. 



This restriction arises because we are working with a relativistic 
theory in which c introduces a fundamental unit of velocity. Non- 
rclativistic self-similar models have the freedom to choose any 
power-law for the radial profile of the angular velocity, including 
in particular a Keplerian profile, but the price for this is that they 
ignore all relativistic effects. 



2.3 The Alfven Critical Surface 

The highest derivative T" in equation 1261 is multiplied by 
a factor of the form 



1 - M T 



-2/i 



= 1- 



n 2 R 2 



(30) 



and so the differential equation is clearly singular when 
this factor goes to zero. The singularity corresponds to the 
Alfven critical surface, also called the light cylinder in some 
of the earlier literature. It is located at the radius 



Ra -w 



(31) 



Note that the Alfven surface occurs at a specific value of 
u — u c and has the shape of a cone. This is in contrast to the 
pulsar magnetosphere problem where it occurs at a single 
radius, thereby motivating the terminology "light cylinder" 
in that case. 

When it = u c the factor given in eq llMOl goes to zero, 
and therefore in order for the solution to be well-behaved the 
rest of the terms in 1261 that do not involve T" should also 
add up to zero. We thus obtain the following two regularity 
conditions at it c : 



T(u c ) 
T'(u c ) 



AT, 
-vM v 



(v-\)s 2 



ttS + l 



1/2 



(32) 
(33) 



The negative sign on the second condition is required for a 
physically sensible solution. 

For a given solution T(u), it is easy to determine 
whether any local stretch of the solution is inside or out- 
side the Alfven surface ; here and elsewhere, we use the 
term "inside" for u < u c and "outside" for u > u c . If 
T(u) > T(u c ) = M" , then we are inside the Alfven sur- 
face, and if T(u) < M v we are outside. Note that, for v < 1, 
the condition 1331 can be satisfied for any value of s. How- 
ever, for v > 1, a solution is possible only if s 2 ^ l/(y — 1). If 
|s| is larger than this limit, the solution cannot pass through 
the Alfven surface. 



2.4 Flow Velocity 

While the azimuthal velocity of a field line passing through 
(R, z) is given by equation 1221 , this does not determine the 
fluid velocity. By assumption, the plasma has no inertia, so 
we are allowed to add any velocity parallel (or antiparallel) 
to the field line without changing the solution in any way. 
There is thus considerable ambiguity as to what we mean by 
the fluid velocity. We could, however, ask what the minimum 
fluid velocity is, i.e., we could choose the parallel velocity 
such that the net velocity of the plasma is a minimum. This 
minimum velocity is the "particle drift velocity," 



1 51 



B 2 /4tt' 



(34) 



where 



S = —E x B 
4ir 



^ [-B$B R R+{B 2 R +Bl 



B^Bzzl 



(35) 



© 2006 RAS, MNRAS 000. HH^ 



Self-Similar Force-Free Wind From an Accretion Disk 5 



is the Poynting flux, E is the electric field given in equation 
Q, and B is the total magnetic field strength. We then have 



c B .Rlc B 



B_Bp_ _ R B p 
' R Cp B ■ 



(36) 



where B p is the magnitude of the poloidal component of the 
magnetic field. 

Another relevant speed is the "energy flow speed," 
which is given by 



v E 



\S\ 
U ' 



(37) 



where U is the electromagnetic energy density given by 



U 



5 [I 



El 



Bi (Q, 2 R 2 



The energy flow speed is related to the particle drift speed 
by 

^ = IT<^' (39) 

which shows that ve is always ^ c. Intuitively, one expects 
a physically sensible solution to have ve , Wmin — > c at large 
distance from the disk, which is equivalent to the condition 
that the fast critical surface is at infinity. This requirement 
is met when the following two conditions are satisfied: 



Btf, 
B„ 



MR 



(40) 



The second condition says that we have to be infinitely 
far outside the Alfven surface, i.e., we should focus on the 
"paraboloidal" solutions discussed in § 3. 

So long as we are inside the Alfven surface, we are guar- 
anteed that D m i n < c. However, once a field line goes outside 
the Alfven surface (QR > c), there is nothing in the problem 
to prevent tWn from exceeding c. This highlights a major 
weakness of the force-free model. Even though the model is 
fully relativistic and causal, it has no way of enforcing v < c 
for the fluid it is meant to describe because it does not ex- 
plicitly make use of the inertia. Note that in any region that 
has « m i„ > c, the electric field strength exceeds the magnetic 
field strength, allowing unlimited electrostatic acceleration 
of charged particles. 

In this paper, we focus our attention on solutions that 
have Umin c for all u of interest, i.e., all R and z spanned 
by a given problem. We call such solutions physically viable 
and consider solutions that violate this condition as "un- 
physical." The restriction to physical solutions plays the role 
of a boundary condition on the solution (see § 12. 51 . From 
equation 13611 we see that requiring the solution to be physi- 
cal is equivalent to requiring a minimum amount of toroidal 
field, or equivalently, a minimum level of enclosed j z (see 
the comment below eg. lilt . 

We should note, however, that even solutions we con- 
sider unphysical under the above strict criterion may have a 
role under some circumstances. For instance, an MHD flow 
with plasma inertia may have a force-free zone that matches 
the corresponding segment of one of our unphysical solutions 
but may then deviate into a non- force- free flow in the region 
where the unphysical solution has « m in > c (e.g., Contopou- 
los 1995a) . A comparison of the force- free solutions discussed 
in this paper with self-similar MHD solutions would thus be 
very instructive. 



2.5 Boundary Conditions 

Let us count the number of degrees of freedom in the prob- 
lem and the boundary conditions we need. Since equation 
1261 is a second order differential equation, it requires two 
boundary conditions. The constants v and M are clearly de- 
termined by the properties of the disk and should be treated 
as externally supplied parameters of the problem. However, 
the other two constants, s and u c , are not free parameters. 
They have to be determined self-consistently in the process 
of solving the problem, i.e., they are eigenvalues. 

Consider first a problem in which field lines cross the 
Alfven surface. To determine the two eigenvalues s and u c we 
have the two regularity conditions 13211 and 1331 . Of the two 
boundary conditions needed by the differential equation, we 
have already discussed one, viz., the trivial condition 11811 . 
The second boundary condition must be set at infinity. For 
this purpose we make use of the velocity constraint v m i n ^ c 
discussed in § 12.41 except that we now require v to be exactly 
equal to c at infinity: 



V C Ju^oo V Rip B J 



= 1. 



(41) 



This condition is equivalent to the statement that we pick 
the physically allowed solution with the minimum level of 
toroidal field. In the language of Michel (1969), it is the 
"minimum torque" (or equivalently, minimum energy) solu- 
tion. Such solutions have the minimum angular momentum 
and energy flux that connect to asymptotic infinity (Kennel 
et al. 1983). 

The logic for picking this particular solution will become 
clearer in § 3.3.3 below. The condition Q410 states that the 
"fast critical surface" or the "light surface," i.e., the point 
at which the fluid speed is equal to c, is located at infinity, 
and not "beyond infinity" (see Okamoto 1974, who calls this 
surface the Alfven point). If we are interested in a system 
that extends to a finite maximum value of u — it max rather 
than to infinity, we may wish to use the same condition 
1411 but applied at u m ax- This would be equivalent to an 
outflowing boundary condition on the force-free flow. 

When we consider a problem in which field lines are 
entirely inside the Alfven surface, we need two boundary 
conditions plus a third condition to determine the eigenvalue 
s. In this case, we have only one boundary condition fea. H8|l . 
and the other two conditions must be set at infinity. Given 
a particular solution such as those discussed in § 3.5, it is 
straightforward to come up with the necessary boundary 
conditions. 



3 ANALYTICAL RESULTS 

3.1 Blandford's Paraboloidal Solution 

When v = 1, the differential equation (12(il has an analytical 
solution: 



T{u) = {l + u) 



2x1/2 



P(R, Z ) 



{R 2 + z 



2x1/2 



(42) 



This is the self-similar version of the paraboloidal solution 
of Blandford (1976). The location u c of the Alfven surface 
is obtained by solving the condition 132L which gives 
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(l-M 2 ) 

The components of the magnetic field are given by 
B R (R, z) = A 



R R(R 2 + z 2 y/ 2 ' 



B,(R,z) = 
B z (R,z) = 



(R 2 + Z 2 ) 1 / 2 ' 



(43) 

(44) 

(45) 
(46) 



The poloidal projection of a field line with footpoint at Rf p 
has the shape 



Z = Jr- (R 2 - R'f 



2R 



fp 



tp) , 



(47) 



which shows explicitly the parabolic geometry of the field 
(hence the name of this solution). 

Note that the poloidal field structure is independent 
of the rotation parameter M as well as the eigenvalue s. In 
order to determine s, we must apply the boundary condition 
14111 . Substituting equations 1441 1461 in equation 1361 1. we 
find 

2 (48) 



«min ■ 

C / u — oo M 



Thus the boundary condition 11 II gives 



s = 2. 



(49) 



Only for this value of s does the fluid velocity u m i n asymp- 
totically tend to c at infinity, (s = —2 is also a solution, 
but this corresponds to a Poynting flux pointed towards the 
disk and is not of interest.) For s < 2, u m i n exceeds c be- 
yond a certain value of u, and the solutions are unphysical. 
For s > 2, w m in is less than or equal to c as u — * oo and it 
would appear that all the solutions with s > 2 are physical. 
Within the framework of the pure self-similar problem, there 
is nothing to indicate that any one of these solutions is more 
physical than the others. However, the more complete anal- 
ysis given by Blandford (1976), which includes a regularity 
condition as R — > 0, shows that only s = 2 is consistent with 
a non-singular solution on the axis. All other values of s that 
appear to be physically valid within the self-similar analysis 
are, in fact, unphysical since they require a singular current 
along the axis. This point will become clearer in § 3.3.3. 

The fact that Br and B z of the paraboloidal solution 
are independent of M means that the poloidal structure is 
independent of rotation. Thus, the presence of the azimuthal 
field B^, however strong it might be, has no collimating or 
decollimating effect. Roughly, one could say that the colli- 
mating hoop stress associated with the toroidal field is ex- 
actly canceled by the decollimating effect of the pressure 
gradient associated with the same field (see Ostriker 1997). 

For the solution with s — 2, we can calculate the 
Lorentz factor 7 m i n corresponding to the minimum fluid ve- 
locity Vmin- Assuming u 3> 1, we find as a function of dis- 
tance from the disk plane 



2M 



(l-M 2 )i/2 \R tp 



1/2 



V2 



R 2 



R 2 p 



R lc 



1/2 



,(50) 



where Rf p is the radius of the foot-point of the field line 
under consideration. Beskin & Nokhrina (2006) derive and 



discuss the asymptotic scaling, ""/mm « R/Ra, for a relativis- 
tic MHD outflow. 



3.2 Michel's Monopole Solution 

When v — 0, the differential equation 126H again has an 
analytical solution: 



T(u) = 1 
P(R,z) = 1 



(l + «2)i/a' 



(R 2 + z 2 y/ 2 ' 

The components of the magnetic field are given by 
R 

B ^ Z ) = (i? 2 +z2) 3/2- 

B^R,z) = --- '* 



(51) 



B z (R,z) = 



c (R 2 + z 2 ] 

z 

(R 2 + z 2 )3/ 2 > 



(52) 



where, as in the paraboloidal case, we have used the bound- 
ary condition 14111 to determine the coefficient of B^. This is 
the monopole solution of Michel (1973a), with the angular 
velocity parameter fl replacing M. 

The monopole solution is not really a disk wind solution 
but rather corresponds to a wind emitted from a rotating 
star. However, by stitching together copies of the monopole 
solution with opposite signs in the two hemispheres it is easy 
to generate a split monopole solution. This has a current 
sheet in the equatorial plane, just like the other solutions 
described in this paper, and is more like a disk wind. In 
both the pure monopole solution and the split monopole 
solution, the poloidal field is purely radial and its geometry 
is independent of the angular velocity of the star Q. Thus, 
as in the paraboloidal solution, rotation induces field lines 
to sweep back (see the expression for B^), but it has no 
tendency to collimate the poloidal field. 

The Alfven surface is at a fixed radius _Ra = O/c and 
takes the form of a cylinder (the "light cylinder"). This is 
not the case for any other value of v. 



3.3 Asymptotic Power-Law Solution at Large u 
for General v 

The two analytic solutions discussed above both have T(u) 
decaying as a power-law at large u; Blandford's paraboloidal 
solution (y — 1) goes as T(u) oc 1/u and Michel's monopole 
solution (y = 0) goes as T(u) oc 1/u 2 . We now consider a 
general value of v and seek a power-law solution of the form 
T(u) oc l/u M with positive \i. For convenience, we will refer 
to such solutions as "paraboloidal" since the field lines have 
the shape of generalized parabolae (see eg. I55H . 

In the limit of large u, we are well outside the Alfven 
surface (recall that "outside" means u > u c ), and the terms 
proportional to M 2 in equation 1261 dominate over the other 
terms. We may thus neglect the subdominant terms and we 
may also replace u 2 + 1 by u 2 . We are then left with the 
simpler equation 

u 2 TT" -u 2 T' 2 /v+(3-2v)uTT' + v(v-l)(l-s 2 )T 2 = 0, (53) 



which has a power-law solution of the form 



T(u) oc u 



1). 



i/^O, 1. 



(54) 



© 2006 RAS, MNRAS OOO.HH^ 



Self-Similar Force-Free Wind From an Accretion Disk 7 



The above solution for /i is not valid when v = or 1 because 
the term with s 2 in equation 1531 vanishes for these cases. 
However, we have fully analytic solutions for precisely these 
values of v, as discussed in the previous subsections. Here 
we are interested in other general values of v. In order to 
be outside the Alfven surface, we require T(u) < M v (as 
discussed earlier), and this is asymptotically possible only if 
[i, is positive. Thus we require s > 1. Using equation 12011 wc 
find that the poloidal shape of the field lines is given by 



z oc R 



s/(s-l) 



R oc z 



(s-l)/s 



(55) 



The shape is vaguely paraboloidal in the sense that z goes 
as a power of R, though the exponent is not exactly 2. 

Given the above solution for fj,, the magnetic field com- 
ponents can be calculated via equations 127H291 and the 
fluid velocity i> m in can be obtained from equation 1361 . It 
is easily verified that any asymptotic power-law solution of 
the form 1541 automatically satisfies the boundary condition 

Vmin = C at U = OO. 

Appendix A discusses the properties of these power- 
law solutions in some detail. Equation 1A141 gives a more 
accurate estimate of the shape of the field line, including 
a first order correction term. Expressions are also given for 
the scalings of the three components of the magnetic field 
and the minimum Lorentz factor 7 m i n . The results for the 
latter generalize the work of Beskin & Nokhrina (2006) who 
consider the perfect paraboloidal case v = 1; however, they 
analyze the MHD problem whereas we consider the simpler 
force-free case. 

To conclude this subsection, we write down the scalings 
of various quantities corresponding to two regimes: v > 1 
and v < 1. These are the key results of this paper. The scal- 
ings are derived in the Appendix and are discussed further 
in § 4. We also provide additional insight into the boundary 
condition at infinity. 



3. 3. 1 Scalings for v ^ 1 

As shown in the Appendix, the case v > 1 is simple. For 
any given v, there is only one value of s for which there is a 
paraboloidal solution; this value is almost exactly equal to 
2/v. The asymptotic form of the solution is 

T(u) « u~ {2 - v \ (56) 
and the poloidal shape of a field line is given by 



(57) 



Rip ~ \R ip ) 

We give only the leading order terms here, and the reader 
is referred to Appendix A. 2 for the next order terms. The 
three components of the magnetic field scale as follows along 
a field line: 

/ „ \ / D \ (4-*/)/2 



\RtpJ 

Mi) 



-M 



R ip 



(2-iz) /2 



Rip 



(3-1/) 



(58) 
(59) 
(60) 



M 



Rip 



(2-f)/2 



M 



R 



R 



Rip Ra 



(61) 



The final expression, R/Ra, is identical to the result ob- 
tained by Beskin & Nokhrina (2006). 



3.3.2 Scalings for v < 1 

When v < 1, paraboloidal solutions are present for all values 
of s > 1, and there is no longer a unique solution. However, 
the "correct" solution is still very close to s = 2/v (§§ 4, 5), 
so the scalings given above for T(u), z/Rt p , Br, B$, B z , are 
all approximately valid. The minimum Lorentz factor now 
has a more complicated dependence. We find 

/ \ (2-i/)/2 p 

- M {ltp) itp <uo > (62) 



7min 



Rii 



v/1 



R 



> wo, 



fp 



(63) 



where u = M~ 1/{1 ~ v) . We note that the Beskin & Nokh- 
rina (2006) scaling 7 m i n ~ R/Ra is valid for small values of 
z, but the scaling is different at large distance. The latter 
regime may be interpreted as the acceleration produced by 
field line curvature, whereas the former regime represents 
acceleration as the result of a "slingshot" effect (V. Beskin, 
private communication) . 



3.3.3 Boundary Condition at Infinity as a Regularity 
Condition on the Axis 

The paraboloidal solutions described above asymptotically 
have £IR 2> c, and their poloidal streamlines are almost pre- 
cisely cylindrical. Let us, therefore, consider a purely cylin- 
drical force-free problem in which Br = 0, and B$, B z and 
Q. are functions only of R. The equilibrium condition (Grad- 
Shafranov equation) for such a flow takes the simple form, 



dBl 



dR R 2 dR 



44 (r^bi 



Asymptotically far from the disk our self-similar solutions 
have i?0 3> B z , and so the first term is negligibly small; in 
fact, for the particular self-similar solution with s = 2/v dis- 
cussed in § 3.3.1, B z is independent of R at a fixed asymp- 
totic 2, so this term vanishes exactly. The remaining two 
terms then give 



B 



^TbUc, 



(65) 



The minimum Lorentz factor varies with z as 



where C is an arbitrary integration constant. For a given 
B z , we have an infinite number of solutions, each with a 
different value of C. This is equivalent to the freedom we 
had in choosing s in §§3.1, 3.2. However, as we now show, 
only one of these solutions is well-behaved. 

Although we have focused on self-similar solutions in 
this paper, we imagine that the solutions are modified inside 
a "core" region at small radii so as to be analytic on the axis. 
The rotation profile, for instance, might behave as a power 
law at large R but we expect it to asymptote to a finite 
constant value as R — * 0. Equation 1651 then shows that 
B^, — * C 1//2 on the axis. However, any solution that has a 
finite Btj, as R — * will have a singular current along the 
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axis, which might be considered undesirable. Therefore, let 
us impose the additional condition that B$ — ► as R — ► 0. 
This can be satisfied only if C — 0, and so we obtain the 
unique solution 



OR 

c 



(66) 



By I65H this relation must be true not only near the axis but 
at all R, including the self-similar region outside the core. 
Thus, analyticity on the axis picks out a unique solution even 
in the far-field self-similar zone. Blandford (1976) reached 
the same conclusion for the specific case v = 1. 

In the self-similar region, our solutions have B$ 3> B z , 
so the quantity on the left-hand side of 16611 is equal to B/B p , 
while that on the right is equal to R/Ra = MR/Rf p . Thus, 
the condition 1661 is the same as 



Rfp B 



(67) 



i.e., the same as the condition I4H . In other words, the 
boundary condition n m i n — > c at infinite distance, which we 
previously argued is a physically motivated boundary condi- 
tion, is in some sense equivalent to the requirement that the 
solution be analytic on the axis. This provides additional 
insight on the boundary condition 1411 . 



3.4 Power-Law Solution Inside the Alfven Surface 

Consider next the region well inside the Alfven surface. Here 
we may ignore the terms proportional to M 2 in equation 
12611 . Let us also assume u 2 3> 1. The differential equation 
then simplifies to 



uT" + (3 - 2v)uT' + v(v - 2)T = 0. 
This has a power-law solution 
T(u) oc tt -M , j-i = —v, 2 — v. 



(68) 



(69) 



Since we assume v > 0, the solution \i = — v corresponds to 
T(u) increasing with u. Asymptotically, this gives a solution 
with a poloidal streamline of the form 



2 oc R° — constant. 



(70) 



Thus, as u — ► oo, the field line asymptotes to a constant 
value of z and the radius R tends to zero. Thus, field lines 
converge radially onto the z-axis. This is clearly unphysical 
since it implies a non-zero magnetic monopole density on 
the axis. Interestingly, the self-similar solutions described 
by Blandford & Payne (1982) are of this kind and are un- 
physical. 

The second solution in 16911 . n — 2—v, has T(u) decreas- 
ing with increasing u (so long as v < 2). This asymptotic 
solution is valid only over the restricted range of u satisfying 
u c ^> u^> 1. For the specific case of a non-rotating solution, 
u c — > oo, and the solution is valid for all u 3> 1. The poloidal 
shape of the field line is given by 



z oc R 



2/(2-") 



(71) 



This corresponds to a perfect parabola when v = 1 (Bland- 
ford's solution), and is a generalized parabola for other val- 
ues of v. 



3.5 Asymptotic Cylindrical Solution 

Consider next solutions in which field lines asymptotically 
become cylindrical with a finite value of R. That is, as 
2 — > oo, the radius R of each field line tends to a finite value. 
Equivalently, as u —* oo, we have both T" and T' tending 
to and T tending to a finite value. Note that these con- 
ditions are different from the quasi-cylindrical asymptotics 
we discussed in § 3.3.3. There, R/Rt p of a field line went 
to infinity as u — > oo, whereas here we consider solutions in 
which R/Rtp remains finite. 

Given the above conditions, most of the terms in equa- 
tion 126H vanish and we immediately obtain 



T(u) 



(y-2) 



,/2 



(72) 



The asymptotic radius Roc of the field line is related to the 
radius of the footpoint by 



Roo 

Rtp~ 



(v-2) 



-1/2 



(73) 



and the asymptotic components of the magnetic field are 
given by 



Br = 0, 

= - 

B z = R v ~ 2 v 



(u-l)(l-s 2 ) 



{v-2 




"l)(l-s 2 )l 




v-2) \ 





1/2 



The azimuthal velocity of the field line is 

-1/2 



"</> x 
C / z-,c 



(v-2) 



(74) 



(75) 



The above expressions give real numbers only if the 
quantity in the square parentheses is positive. If v > 1, 
this requires s 2 > 1, and if v < 1, it requires s 2 < 1. These 
are necessary conditions for the existence of cylindrical so- 
lutions. 

Note that asymptotically cylindrical solutions may be 
either inside or outside the Alfven surface. Whether a par- 
ticular solution is inside or outside is determined by whether 
-Roo is less than or greater than _Ra, or equivalently, by 
whether v$/cis less than or greater than unity. If the quan- 
tity in square parentheses is less than 1, then the solution 
is outside the Alfven surface. Clearly, this is always the case 
when v < 1. When v > 1, there is no apparent restriction, 
though in practice only solutions inside the Alfven surface 
are present (see § 4). 



3.6 Conical Solution 

These are solutions in which T goes to zero at a finite value 
of u — Uf, so that field lines asymptotically have a conical 
shape with z/R = Uf as z, R — > oo. The solution does not 
exist for R < z/uf, i.e., inside the cone. Thus, for force bal- 
ance, one must supply the required pressure at the surface 
of the cone. Because this is somewhat unphysical, and also 
since w m i n > c as u — > uf, we do not consider these solu- 
tions interesting. Nevertheless, we discuss them briefly for 
completeness. 

As stated above, we have T — * as b - * Uf. In this 
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limit, the dominant terms in equation 1261 are the first and 
second terms inside the square brackets. Focusing on these, 
the solution must satisfy 

vTT" = T' 2 , (76) 
which has a power-law solution of the form 
T{u) oc (u f -u) u/(v - 1} , v>\. (77) 
This is the conical solution, which is present only for v > 1 . 

4 NUMERICAL SOLUTIONS 

Before describing the numerical solutions, let us recall the 
properties of Blandford's (1976) paraboloidal solution (y — 
1, § 3.1). 

(i) The poloidal structure of field lines is perfectly 
parabolic. 

(ii) Rotation has neither a collimating nor a decollimating 
effect on the poloidal configuration of field lines (see Fig. 1). 

(iii) For values of s 2, the solutions are physically ac- 
ceptable in the limited sense that n m j n ^ c at all points (solid 
and dashed lines in Fig. 2). For s < 2, the solutions have 
«min exceeding c beyond a certain value of u. (Technically, 
we should replace s by |s| in these statements, but we are 
interested only in solutions with positive s.) 

(iv) Among solutions with s 2, the particular solution 
with s = 2 has n m j n asymptotically tending exactly to c 
(fast critical surface at infinity), and it is also the only so- 
lution that can be matched to an analytic core on the axis 
(§ 3.3.3). In addition, it is the solution with the minimum 
torque among all solutions that have w m i n $C c. Thus, this 
unique solution is the "correct" solution to the self-similar 
problem. In a sense, the constraint on the fast surface and 
the analyticity argument provide a physical motivation for 
Michel's (1969) minimum torque condition. 

(v) For the solution with s = 2, the Lorentz fac- 
tor increases with distance along the z-axis as 7 m i n ~ 
M(z/R tp ) 1/2 feci. 1501 Beskin & Nokhrina 2006). 

As we show below (and also in the Appendix), the 
paraboloidal solution v — 1 acts as a watershed in solution 
space — solutions on one side (y > 1) have very different 
properties compared to those on the other side (y < 1). 

4.1 Solutions with v > 1 

We have numerically scanned the solution space of equation 
12611 for the particular choice of parameters v — 1.25 and 
M — 0.1. Keeping these parameters fixed, we solved for T(u) 
corresponding to a range of values of the eigenvalue s. Then, 
using the numerical solution we calculated the components 
of the magnetic field, the minimum fluid velocity, the mini- 
mum Lorentz factor, etc. We identified which solutions are 
physically acceptable, i.e., have v m i n ^ c for all u (techni- 
cally, out to the maximum u to which we integrate, which is 
~ 10 12 ), and which are not, which solutions have power-law, 
i.e., generalized parabolic, asymptotics (eq. 1551 and which 
have cylindrical asymptotics as u — * oo, and we also flagged 
any special solutions that had n m j n — > c as u — > oo. 

As per the discussion in § 2.3, we first checked whether 
s 2 ^ l/(y — 1). If so, we solved for u c , i.e., the position of 




12 3 4 

log (1+R/RJ 



Figure 1. Shows the poloidal structure of field lines for the 
paraboloidal solution, v = 1. The calculations correspond to 
M = 0.1 and eight values of s: 0.5, 1.0 ■■■ 4.0. All the solu- 
tions have exactly the same poloidal structure (§ 3.1). Contrast 
this with Figs. 3 and 5. 




2 4 6 8 10 

log(l+z/R, p ) 




2 4 6 8 10 



log(l+z/R, p ) 

Figure 2. Shows the variation of the minimum Lorentz factor 
7min and the minimum velocity /3 m in — ^min 

/c as a function 

of distance z from the disk plane for the paraboloidal (v = 1) 
solutions of Fig. 1. The thick solid line is the physically interesting 
solution with s = 2, whose velocity asymptotically approaches c 
at infinity. The dotted lines correspond to solutions with s < 2, 
which attain minimum velocities exceeding c. The dashed lines 
correspond to solutions with s > 2, whose velocities remain below 
c at all z. 
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Figure 3. Shows the poloidal structure of field lines for self- 
similar solutions with v = 1.25. The calculations correspond to 
M = 0.1 and s = 0.5, 1.0, 1.5, 1.59, 1.599, 1.5999, 1.59996 (dotted 
lines), s = 1.5999635948 (heavy solid line), s = 1.7, 1.8, 1.9, 2.5, 
3.0, 3.5, 4.0 (dashed lines). The heavy solid line is the physically 
interesting solution that satisfies the boundary condition 1411 . 



the Alfven surface, by integrating equation 1261 and apply- 
ing the conditions 1181 . 13211 . 13311 . In practice, we found it 
convenient to assume a value of u c and integrate the differ- 
ential equation 1261 from the Alfven surface u = u c , where 
132H and 133H provide initial conditions, back to u — 0. At 
u — we checked whether the solution satisfied the condi- 
tion T(0) = 1. If not, we tried a different value of u c and 
iterated until the solution converged. We then integrated 
equation 1261 outside the Alfven surface to large values of 
ii > ii c to complete the solution. 

If s 2 > 1/0-1), we looked for a solution that lives en- 
tirely inside the Alfven surface. The only physically accept- 
able solutions of this kind are those with cylindrical asymp- 
totics. Therefore, we used the discussion in § 3.5 to assign 
the necessary boundary conditions at infinity. By construc- 
tion, all of these solutions have v m in < c. 

Figure 3 shows the poloidal structure of the field lines 
for various solutions, and Figure 4 shows the corresponding 
variations of /3 m i n = «W n /c and 7 m i n with distance. The 
thick solid lines correspond to the unique solution that has 
Wmin — » c as u — » oo. This solution corresponds to s « 1.6 = 
2/v (see Appendix AT). The poloidal shape of field lines is 
asymptotically of the form z oc R 2 ^ 2 ~ v ^ (see eg. 15511 . i.e., a 
generalized parabola. 

Solutions with s on either side of the above critical value 
are "unstable" (Appendix AT; see also Contopoulos 1995a). 
When s <2/v (shown by dotted lines in Figs. 3 and 4), the 
solutions have fluid velocities exceeding c and are unphysi- 
cal. They are all asymptotically conical (see § 3.5). Solutions 
with s > l/(v — l) 1 / 2 are entirely inside the Alfven surface 
and have cylindrical asymptotics (dashed lines in Figs. 3, 4). 
These solutions all have n m j n < c as u — * oo. Over the range 



Figure 4. Shows the variation of the minimum Lorentz factor 
7 m i n and the minimum velocity /3 m i n = f m i n /c as a function of 
distance z from the disk plane for the v = 1.25 solutions of Fig. 
3. 

2/v < s < l/(v — 1) 1/2 , field lines go through the Alfven 
surface, move out to a maximum radius and then turn back 
towards the Alfven surface. The solutions become singular 
when they attempt to recross the Alfven surface. 

4.2 Solutions with v < 1 

Figures 5, 6 show results corresponding to the parameter 
choice of v — 0.75 and M = 0.1. As before, we have calcu- 
lated solutions for a range of values of s and cataloged their 
properties. 

For s ^ 1 (the first five dotted lines in Figs. 5, 6), there 
are no asymptotic power-law solutions available (see the dis- 
cussion in § 3.3). These solutions are asymptotically cylin- 
drical, superposed with a decaying oscillation, and they are 
unphysical (y m i n > c). For all values of s > 1, we have 
power-law paraboloidal solutions, and these are "stable" in 
the sense described in Appendix AT. 

Over the range 1 < s < 2.8146, the solutions have 
^min > c somewhere within the range of integration over 
u and are thus unphysical. The thick solid lines in Figures 
5, 6 correspond to s — 2.8146, the lowest value of this pa- 
rameter for which we can find a physically allowed solution 
extending over a large range of u (our integrations extended 
out to a maximum u ~ 10 12 ). As per equation 155H . the 
poloidal shape of field lines follows z oc R 155 . We note that 
this critical solution has nearly the same scalings as for the 
v = 1.25 solution described in the previous subsection, viz., 
s w 2/v = 2.67, z oc _R 2 /< 2 -") = R 1(i , but it is slightly dif- 
ferent. The difference arises because of the finite value of M 
that we have selected. In the limit of small M, it turns out 
that the critical solution satisfies exactly the same scalings 
as the v > 1 solutions described in the previous subsection 
(see § 4.3 below). 
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Figure 5. Shows the poloidal structure of field lines for self- 
similar solutions with v = 0.75. The calculations correspond to 
M = 0.1 and s = 0.2, 0.4 •■ • 2.8 (dotted lines), s = 2.8146 (heavy 
solid line), s = 2.9, 3.0, 3.2, 3.4, 4.0, 5.0, 10.0 (dashed lines). 
The heavy solid line is the physically interesting solution. It is 
the solution with the smallest value of s, i.e., the least amount 
of toroidal field (minimum torque), that satisfies the boundary 
condition 1411 . 




2 4 6 8 10 



log(l+z/R, p ) 

Figure 6. Shows the variation of the minimum Lorentz factor 
7min and the minimum velocity /3 m i n = ^min /c as a function of 
distance z from the disk plane for the v = 0.75 solutions of Fig. 
5. 




0.5 1 1.5 2 



v 



Figure 7. Demarcates the different solution regimes for self- 
similar force-free winds from a thin disk. The regimes are indi- 
cated as a function of the similarity index v and the field sweep- 
back parameter s. The results correspond to the limit M < 1. 
The thick line is the locus of solutions with s = 2/u. These are 
the physically interesting minimum torque solutions; at each u, 
this is the lowest value of s for which a physically valid solution 
with v m i n < c is possible. The various regimes demarcated by 
dashed lines are described in § 4.3. 

For s > 2.8146 and up to a second critical value ~ 3.8 
(the first four dashed lines in Figs. 5, 6), all the solutions 
are physically acceptable, and surprisingly all of them have 
Vmin — » c as u — > oo. Thus, there is a continuous family of 
solutions all of which satisfy the boundary condition 14111 . 
However, all these solutions have larger toroidal fields, i.e., 
larger torques, than the particular solution with s = 2.8146 
described in the previous paragraph, and their Lorentz fac- 
tors increase more slowly with increasing z. Finally, for 
s > 3.8 (the last three dashed lines), the solutions revert 
to being unphysical and the fluid velocity exceeds c over 
some range of it. 

We should note that the results of this subsection de- 
viate from those of Contopoulos (1995a) who states that 
there are no physically viable solutions for v < 1 (see his 
discussion of solutions of type 1). 



4.3 Survey of Solution Space 

The discussion in the previous two subsections was for two 
specific values of v and for a single choice of M. We now 
briefly put these results in a more general context. 

We begin by considering the case in which the rotation 
parameter M is very small. Figure 7 shows for this case the 
solution space of self-similar force-free disk winds as a func- 
tion of the stream function index v and the field sweep-back 
parameter s (compare with Fig. 4 in Contopoulos 1995a, 
noting that s — * —Hoc n and v — * x). The thick solid line 
corresponds to the relation s = 2/u. 
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For v > 1, paraboloidal solutions are available only on 
the thick solid line. Below the line, we have unphysical con- 
ical solutions. Above the line, within a narrow triangular 
region, we have no solutions at all. In this region, solutions 
cross the Alfven surface and then try unsuccessfully to re- 
cross the Alfven surface (§4.1). The solutions are thus unable 
to reach infinity. Above this forbidden zone, there is an ex- 
tended region of parameter space where there are asymptot- 
ically cylindrical solutions, all of which are physically consis- 
tent. The main difference between these solutions and those 
on the thick solid line is that the cylindrical solutions have 

l>min < C as U — > OO. 

When v = 1, paraboloidal solutions with physically al- 
lowed velocities are available for all points above and on the 
thick solid line. However, the point s — 2, which is exactly 
on the line, is the only solution that is analytic on the axis 
(§ 3.3.3), and we have argued that it is therefore the correct 
solution. 

For v < 1, the situation is quite interesting. In the limit 
M< 1 that we are considering here, physically acceptable 
paraboloidal solutions are available for all points above and 
on the thick solid line. Moreover, all these solutions have 
Vmin — » c for large u and all are analytic on the axis. How- 
ever, the solutions with s = 2/u (the thick solid line in Fig. 
7) are still special in that these solutions have the mini- 
mum torque among all the physical solutions and also the 
most rapid acceleration outward. Previously, we associated 
the condition v m - m — c at infinity with analyticity on the 
axis and also with the minimum torque condition. Now we 
find a continuum of solutions that have u m i n = c at infinity 
and are analytic, and it is only the requirement of minimum 
torque that picks out a unique solution. Following Michel 
(1969), we believe the minimum torque solution is the phys- 
ically relevant solution. Points below the thick solid line in 
Fig. 7 are unphysical since u m i n exceeds c. These unphysical 
solutions are either paraboloidal or cylindrical, as indicated 
in Figure 7. 

When M is not arbitrarily small, the space of phys- 
ically allowed solutions is virtually unchanged for v ^ 1. 
However, for v < 1, the space of allowed solutions shrinks. 
This is shown by the four dotted lines in the left panel 
of Figure 8, which correspond (from the left) to M — 
0.0001, 0.001, 0.01, 0.1. In each case, only the region to 
the right of the dotted line is physically allowed, while the 
region to the left (even if it is above the thick solid line) 
has Vmin > c. Also, for each value of v, the lowest value of s 
for which a physically allowed solution is present (the thick 
segment of the dotted line) is the one with the minimum 
torque and the most rapid acceleration. In effect, this solu- 
tion takes on the role of the solution on the thick solid line 
when M < 1. 

We should note that the above statements about physi- 
cally allowed and disallowed solutions are based on integrat- 
ing our solutions out to a large value of u ~ 10 12 (beyond 
this, we are not confident of the accuracy of the numeri- 
cal integration.) However, in practice, one rarely requires a 
force-free solution to be physical out to such large distances 
since even very high-a relativistic MHD flows are likely to 
feel the effects of inertia well before this. It is therefore inter- 
esting to look for force- free solutions that are physically rea- 
sonable (f m in < c) out only out to a modest value of u. The 
right panel in Figure 8 shows the allowed region in param- 




Figure 8. Shows further details of self-similar solutions for v < 1. 
[Left] The thick solid line is the same as in Fig. 7. In the limit 
M -C 1, points on and above this line correspond to physically 
consistent solutions and points below are inconsistent. The four 
dotted lines show how the region of consistent solutions shrinks 
for finite values of M; from the left, the lines correspond to 
M = 0.0001, 0.001, 0.01, 0.1. For each choice of M, the re- 
gion to the right of the corresponding dotted line is the region 
of consistent solutions. All solutions were numerically integrated 
from it = to M max ~ 10 12 to determine whether or not they 
are consistent. [Right] The thick solid line is again the same 
as in Fig. 7. The dotted lines show the effect of reducing the 
value of n max for M = 0.1; from the left, the lines correspond 
to u max = 10 3 , 10 5 , 10 7 , 10 9 . For each choice of Umax, the re- 
gion to the right of the corresponding dotted line is the region of 
consistent solutions. 



eter space for different ranges of integration; from the left, 
10 3 , 10 5 , 10 7 , 10 9 . Clearly, solutions are allowed 
over more and more of the parameter space as we reduce 
the integration range. However, we still have the situation 
that, for each v and integration range, the smallest value 
of s for which a physically allowed solution is available (the 
thick segment of the dotted line), i.e., the minimum torque 
solution, is the one with the largest acceleration. 

We now come to a key question. For a given choice of 
v and M, i.e., for given self-similar boundary conditions at 
the equatorial plane, do we know which of the many solu- 
tions described here is the correct solution? By correct, we 
mean the solution that would be picked out by a real sys- 
tem. We have presented a number of arguments to suggest 
that we should pick the solution with t> m in — > c since this 
condition appears to be equivalent to analytic behaviour on 
the axis. But this condition alone is insufficient when v < 1 
since there are many solutions satisfying the condition. In 
the latter case, we have suggested that the correct solution 
is the one with the minimum torque. The numerical simula- 
tions described in §5 are designed specifically to verify these 
assertions. 
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Figure 9. Variations of the three field components, Br, B$, 
B z , and the minimum Lorentz factor 7 m i n as a function of the 
coordinate z along a field line. The index v is equal to 1.25, and 
three solutions are shown, corresponding to M = 0.1, 0.01, 0.001, 
with s = 1.5999635948, 1.59999998, 1.60000, respectively. The 
solid lines show the numerical results and the dotted lines show 
the approximate scalings given in § 3.3.1. 

4.4 Numerical Verification of Scalings 

In §§ 3.3.1 and 3.3.2, we wrote down scaling relations for 
a number of quantities of interest. We have verified these 
scalings by comparing the analytical results with numerical 
solutions. 

Figure 9 shows results for u — 1.25 and three values of 
M: 0.1, 0.01, 0.001. Shown are the three components of the 
magnetic field, Br, B^,, B z , and the minimum Lorentz factor 
7min along a field line. The solid lines correspond to the 
numerical results and the dotted lines show the analytical 
scalings. As discussed in Appendix A. 2, the analytical results 
ignore numerical coefficients of order unity and so there are 
small vertical offsets between the solid and dashed lines. 
Apart from this, the agreement is good. 

Note that the poloidal components of the field, Br and 
B z , are essentially independent of M. Thus, rotation has no 
effect on the poloidal structure, just as in the paraboloidal 
case discussed in § 3.1. The toroidal component of the field, 
B^, does depend on rotation. In fact, its strength is directly 
proportional to the rotation parameter M. 

The manner in which the Lorentz factor scales with z, 
viz., 7 m in oc z^- 2 ^"^ 2 , is independent of M, but the actual 
value of 7 m in at any given z does depend on M. This is be- 
cause, with decreasing disk rotation, the acceleration of the 
wind starts at a progressively larger value of z. This is not 
surprising. Acceleration is effective only outside the Alfven 
surface, because only there does B,p become the dominant 
component of the field. When the disk rotates slowly, the 
Alfven surface is crossed at a larger value of z. 

Figure 10 shows results for v — 0.75 and the same three 
values of M. For each M, we have selected the lowest value 




Figure 10. Similar to Fig. 8, but for v = 0.75. The three 
solutions correspond to M = 0.1, 0.01, 0.001, with s = 
2.8146, 2.67366, 2.66710, respectively. The solid lines show the 
numerical results and the dotted lines show the approximate scal- 
ings given in §3.3.2. 

of s for which we could obtain a physically valid solution. 
As discussed earlier, this is the solution with the minimum 
torque and the largest acceleration. The poloidal structure 
of the field shows some variation with M, but the changes 
are quite small, so we again find that the poloidal structure 
is effectively independent of rotation. 

The variation of the Lorentz factor with z is rather in- 
teresting in this case. We see that there are two regimes: a 
regime of rapid acceleration soon after a field line crosses 
the Alfven surface, and a regime of slower acceleration far- 
ther out. The scalings given in § 3.3.2 for 7 m i n agree quite 
well with the numerical results. Interestingly, the asymptotic 
Lorentz factor at large z is somewhat larger for a slower spin- 
ning disk than for a rapidly spinning disk. This is seen both 
in the numerical results and in the scaling relations. 



5 TIME-DEPENDENT NUMERICAL 
SIMULATIONS 

In this section, we discuss the results of time-dependent nu- 
merical simulations of force-free nearly self-similar winds. 
The genera l relativisti c force- f ree electrodynamics code de- 
scribed in iMcKinnevI d2006aT) is used to evolve the ax- 
isymmetric force-free equations of motion in a flat space- 
time in spherical polar coordinates. This code is an exten- 
sion of the general r elativistic MHD scheme called HARM 
jGammie et all 2003 ) , such that the force-free equations are 
written as a general set of conservation equations. In the 
ideal force- free degenerate limit, the force-free equations can 
be written in a coordinate basis in conservative form as 

= - JLfy^T/) + S(T), (78) 
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where the index i runs over all three spatial dimensions and 
the index j runs over the two poloidal spatial dimensions 
and T is the stress-energy tensor given by 



A a u , 1 j 2 
= o ir u H — b 



ft" 6", 



(79) 



where 6 M is the 4-magnetic field in the frame moving at 
the drift 4-velocity given by u M and the metric g^" corre- 
sponds to Minkowski space-time with a determinant of g. 
The quantity S(T) is a source term that depends on the co- 
ordinate system. These three evolution equations implicitly 
determine the evolution of the electric field. The magnetic 
field is governed by the induction equations given by 



9t^ 



IB') 



d 



(80) 



where B l is the lab-frame magnetic field that obeys the 
solenoidal constraint given by 



dx iV 



0. 



(81) 



For more details see [ Kom issarov |2002b , 
iGammie etafl J2003ft : iMcKinnevI j2006ah . 

The force-free version of HAR M has been success fully 
used to study pulsar magnetospheres jMcKinnevl2006bl) . the 
Blandford-Znajek split-monopole solution (corresponding to 
v = 0, McKinney 2006a), the Blandford-Znajek paraboloidal 
solution (corresponding to v = 1), and a B landford- 
Znajek -type solution with v = 3/4 iMcKinnev fc Naravanl 
2006b). The MHD version of HARM has been success- 
fully used to study accretion flows around rotating black 
holes where nearly force-free jets self-consist ently form and 
are magnetically accelerated to 7 ~ lOjGanrmie^jl^d] 
| 200l 120041 iMcKinney fc Gammid 120041 ; IMcKinnevI 12006c; 
iMcKinnev fc Naravanll2006al) . 



5.1 Modelling the Self-Similar Solution 

The divergence of the field strength near the polar axis and 
at spherical r — in the self-similar models is problem- 
atic for a time- dependent numerical code. Indeed, as we dis- 
cussed in § 3.3.3, we assume that the self-similar solution is 
modified in a core region at small cylindrical radii so as to 
maintain analyticity around the axis. For analogous reasons, 
we screen out the divergent region of the numerical solution 
near spherical r = in the simulations by introducing an 
artificial "star." The star is centered at r — with unit ra- 
dius and is endowed with a constant angular velocity on its 
surface equal to the angular velocity in the self-similar disk 
model at R — 1. The poloidal field at the surface of the star 
is chosen to be the same as that given by the self-similar 
solution on the spherical polar surface r = 1. 

The computational grid is chosen so that at large 
radii the grid follows the collimating field lines to ensure 
good resolution. The grid is composed of an arbitrary two- 
dimensional space with coordinate directions {xi,X2}- The 
xi grid lines are mapped to the spherical polar radius r ac- 
cording to 



r = Ro + e 11 , 

where _Ro = —3, n = 
from Ri n = 1 to R out 



(82) 

10, and the radius of the grid goes 
= 10 4 in arbitrary units of length. The 



X2 grid lines are mapped to the spherical polar 9 according 
to 



tt\ f l +ta,n- 1 [h(x 2 - 1/2)] 
2J \ tan- 1 ^] 



where 



h(r) = 



r - r 
ri 



(83) 



(84) 



and we choose ro = 0, ri = 10. The index a is chosen 
to be the same as the index in the scaling of the opening 
angle of a field line at large radii in the self-similar solution, 
8j oc r~ a . From the analytic asymptotic solution given in 
equation I54H . we find that a — l/s. 

The resolution for all models is chosen to be 256 x 128. 
The solutions are well-converged compared to low resolution 
models except very close to the poles where the coordinate 
singularity of spherical polar coordinates causes minor arti- 
facts that do not affect the results. 



5.2 Obtaining a Force- Free Stationary Solution 

Steady state force-free numerical solutions with no disconti- 
nuities or surface currents above the disk surface are found 
by choosing boundary conditions determin ed by an analysis 
of the Grad- Shafranov equation (see, e.g.. lBogovalovlFl997l : 
iBeskinl Il997fr . For solutions that pass through the Alfven 
surface at some radius, one is required to fix the magnetic 
field component perpendicular to the conductor (B r for the 
star and Be for the disk) and to specify two other constraints 
at the star or disk, viz., E<j, = 0, and the value of Qf, the 
field line angular velocity. 

For axisymmetric, stationary solutions, the frozen-in 
condition of ideal MHD implies that the field line velocity 
Vi is completely determined by the field Bj and field rota - 
tion frequency Qf (see equation 46 in IMcKinnevI l2006af) . 
Thus, during the simulation the 3- velocity at the stellar sur- 
face and on the disk equatorial plane is set to agree with 
this condition. Such a 3-velocity is generally time-like for 
points inside the Alfven surface, but outside this region the 
3-velocity can sometimes be space-like and unphysical (this 
is analogous to v m i n > c which is discussed extensively in 
previous sections). In the event that the 3-velocity becomes 
space-like during the simulation, the Lorentz factor is locally 
co nstrained to a fixed large value as described in section 2.5 
of IMcKinnevI J2006al) . This safety feature is necessary to 
handle the violent evolution seen early in the simulations, 
but it is usually not activated once the solution approaches 
a stationary state. 

In the work described here, the initial conditions in the 
disk are chosen to correspond to the chosen self-similar so- 
lution, described by the parameters v, M and s. The only 
difference is that we initially set B<f, — and allow the code 
to develop whatever toroidal field it wishes. During the time 
evolution, Be at the disk is held fixed at its initial self-similar 
value, and the field angular velocity £If is set according to 
the self-similar model. On the star, the self-similar value of 
B r is held fixed within r = 1 and flp(r = 1) is set equal to 
Qf(R = 1, z = 0). Thus the stellar and disk values of Qf 
match at R = 1, 2 = 0. 

The initial non-rotating state is far from steady state, so 
the model undergoes violent non-stationary evolution. Even- 
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Figure 11. Shows for u = 1, M = 0.25, the run of various quan- 
tities as a function of z following a particular field line with foot- 
point at ijfp = 2.0. The four panels correspond to (i) the opening 
angle of the field line (9j = tan -1 R/z), (ii) the minimum Lorentz 
factor (7 m i n ), (iii) the toroidal field strength (Bj,), and (iv) the 
pitch angle (a p itch = tan - 1 Bp/\BA"). In each panel, the solid 
line corresponds to the final converged solution obtained from 
the time-dependent numerical simulation, the dotted line corre- 
sponds to the analytical self-similar solution with s = 2, and the 
dashed line corresponds to the analytical solution with s = 4. 
It is clear that the numerical solution agrees very well with the 
s = 2 analytical solution. This is the smallest value of s for which 
a physically consistent (f m i n ^ c) solution is available, and it is 
the minimum torque solution. 



tually, however, it relaxes to a steady state. All solutions 
thus found are necessarily stable to Eulerian axisymmetric 
perturbations. 



5.3 Models with v = 1 

The self-similar mod el with v = 1 corresponds to the 
paraboloidal model of iBlandfordl <ll976t) . As we have seen, 
the poloidal structure of the field is independent of M and 
s. We initialize the simulation with the analytical poloidal 
solution and set the angular velocity profile in the disk to 
correspond to the desired value of M, viz., M = 0.25. From 
the discussion in the previous sections, we know that phys- 
ically allowed solutions are available for all s ^ 2, of which 
the solution with s — 2 is the one that (i) has v m in = c 
at large distance, (ii) is analytic on the axis and (iii) has 
minimum torque. However, recall that is set equal to 
initially (corresponding to s — 0) and we allow the code to 
evolve to whatever value of s it chooses. Which value of s 
does the time-dependent code select? 

Figure HTI shows the results of the numerical model com- 
pared against the self-similar analytical model. The plot 
shows the dependence of various quantities along the partic- 
ular field line that starts at a foot point radius of R{ p — 2.0. 
The top left panel shows the field line's angle away from the 
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Figure 12. Contours of the stream function P showing the 
poloidal structure of field lines for u = 1, M = 0.25. The thin 
solid lines show the numerical solution and the dotted lines show 
the analytical self-similar solution for s = 2. Note the good agree- 
ment. The thick lines are explained in the caption to Fig. 1131 



polar axis as a function of distance from the disk plane. 
This quantity is purely a function of the poloidal struc- 
ture of the field lines, and it agrees very well with the 
self-similar solution. The other panels show the minimum 
Lorentz factor, the toroidal field strength, normalized such 
that B Z (R = 1, z = 0) = 1, and the pitch angle of the field 
line defined as 



Ipitch 



tan 



Bp 

\b2 



(85) 



where B p is the poloidal field strength. For comparison, the 
dotted lines show the analytical results for the case s = 2 and 
the dashed lines show the corresponding results for s = 4 (as 
a counter-example). It is clear that the numerical simulation 
converges to the model with s = 2, for which we have 



Qpitch 



tan 



i / 



(86) 



at large radii. 

Figure IT21 shows the simulation results for the poloidal 
field geometry and the location of the Alfven surface. Over- 
lapping the simulation results is the analytical self-similar 
field geometry and the location of the self-similar Alfven 
surface. There is excellent agreement for the Alfven surface 
and reasonable agreement for the disk field lines. Because the 
simulation has a rigidly rotating star at the center, there is 
a second branch of the Alfven surface corresponding to the 
stellar rotation. The position of this branch is in reasonable 
agreement with the analytical estimate, with minor varia- 
tions due to under-resolution of the Alfven surface at large 
radius. 

Figure lTol shows the inner region of Figure [T2l in order to 
show how the components of the Alfven surface associated 
with the star and the disk merge into a single Alfven surface. 
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Figure 13. Close-up of the central regions of Fig. Q21 highlight- 
ing the location of the Alfven surface, shown by the thick solid 
line. Ignoring the disk, the uniformly rotating star would have 
its Alfven surface at R = 4, shown by the thick vertical dotted 
line. The pure analytical self-similar disk wind solution has its 
Alfven surface at a fixed angle of 8 = 0.49 radians away from 
each polar axis, as shown by the thick sloping dotted lines. It is 
seen that the Alfven surface in the numerical solution smoothly 
interpolates between the two analytical Alfven surfaces. 
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Figure 14. Similar to Fig. 1111 but for u = 0.75. The numerical 
model (solid lines) agrees well with the analytical self-similar so- 
lution with s = 2.7 (dotted lines), but not with the self-similar 
solution with s = 5.0 (dashed lines). The former value of s is ap- 
proximately the smallest for which a physically consistent solution 
("min ^ c) is possible, i.e., it is the minimum torque solution. 



This figure also shows that the field lines are in excellent 
agreement at small radii. 



5.4 Models with v = 0.75 

Unlike the case of v — 1, self-similar models with v — 0.75 
have different poloidal field geometries for each value of M 
and s. The values of Be and Q.f in the disk are specified 
by the boundary conditions, i.e., the parameters v and M, 
and are of course independent of s. For this study we choose 
M = 0.1 and initialize the simulation with the poloidal solu- 
tion corresponding to two values of s: 2.8146, corresponding 
to the minimum torque condition, and 5.0, a much larger 
value. Which value of s will the time-dependent numerical 
simulation choose? 

Figure 1141 shows the same information as Figure 1111 
but for v = 0.75. Regardless of which solution we use to 
initialize the simulation, the final state selected by the time- 
dependent code corresponds to a solution with s « 2.7. This 
is close enough to the special value of s = 2.8146 that we 
claim the final solution is the minimum torque solution, in 
agreement with the proposal of Michel (1969). As discussed 
in section 1431 for finite M and for s below a critical value, 
there is a finite radius beyond which the solution becomes 
unphysical. Indeed, s « 2.7 leads to w m in = c at a finite 
radius somewhat smaller than the size of the simulation box. 
The simulation does reach v m i n = c for r ~ 10 3 , although 
the numerical method would have to be improved to verify 
that this is not just numerical error. 



5.5 Models with v = 1.25 

Like self-similar models with v = 0.75, models with v = 1.25 
have different poloidal field geometries for each value of M 
and s, while the values of Bg and Qf in the disk are the same 
for each s. For this study we choose M — 0.4 and initialize 
the simulation with the self-similar poloidal solution corre- 
sponding to two values of s: 1.6, which corresponds to the 
paraboloidal solution on the thick solid line in Figure 7, and 
s — 2.5, which corresponds to an asymptotically cylindrical 
solution. Again, the interesting question is which solution 
will be picked by the time-dependent numerical simulation. 

Figure H^l shows the same information as Figure [TH but 
for v = 1.25. Clearly the time-dependent code has chosen a 
solution close to the paraboloidal solution with s w 1.6. 
Even when we start the simulation with the cylindrical so- 
lution, the field lines spontaneously open out and become 
paraboloidal with s = 1.6. The final solution is the mini- 
mum torque solution for this problem. 

Note that the behaviour at large radius of the time- 
dependent solution that starts with s — 2.5 is difficult to fol- 
low because the solution undergoes quite a large change from 
cylindrical to paraboloidal streamlines. Thus one under- 
resolves either the initial solution or the evolved final solu- 
tion at large radii. For such models we performed two sim- 
ulations with a grid that resolves one or the other type of 
field geometry, and found that both resulted in the same 
final solution with s ~ 1.6. 
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Figure 15. Similar to Fig. 1111 but for v = 1.25. The numerical 
model (solid lines) agrees well with the analytical self-similar so- 
lution with s = 1.6 (dotted lines), but not with the self-similar 
solution with s = 2.5 (dashed lines). The former value of s is the 
smallest for which a physically consistent solution (v m i n $C c) is 
possible, i.e., it is the minimum torque solution. 

5.6 Summary of Time-Dependent Simulation 
Results 

In summary, the time-dependent force-free numerical simu- 
lations reveal the following: 

(i) For a given choice of v and M, a self-similar force-free 
system spontaneously seeks the smallest value of the eigen- 
value s for which a physically consistent solution (v m i n ^ c) 
is possible throughout the finite numerical grid. Since s mea- 
sures the strength of the toroidal field, an equivalent state- 
ment is that the system chooses the configuration that has 
the least sweepback of field lines in the toroidal direction, or 
the minimum torque. 

(ii) The final converged numerical solution matches very 
well the analytical self-similar solution that corresponds to 
this value of s. There is close agreement in the profiles of 
the components of the magnetic field, Br, B^, B z , and the 
minimum Lorentz factor 7mm- 



6 SUMMARY AND DISCUSSION 

With a view ultimately to understanding the nature of MHD 
outflows from relativistic accretion disks around spinning 
black holes, we have considered in this paper the much sim- 
pler problem of a self-similar force-free outflow from an in- 
finitely thin rotating equatorial disk. The problem is mathe- 
matically described in terms of two dimensionless numbers: 
(i) v, the radial index of the stream function P on the equa- 
torial plane (eq. 1141 . which is defined such that the mag- 
netic field scales as R v ~ 2 on this plane; (ii) M, a rotation 
parameter, defined such that the azimuthal rotation velocity 
on the equatorial plane is Mc, independent of R (eg. I21H . 



Given these two parameters, the self-similar problem nat- 
urally leads to two dimensionless eigenvalues: (i) s, which 
describes the degree of sweepback of the magnetic field lines 
at the equatorial plane (eq. 12311 : (ii) u c , which gives the value 
of z/R at the Alfven surface (also called the light cylinder 
in the pulsar magnetosphere problem). The eigenvalues are 
obtained by solving the differential equation 1261 with ap- 
propriate boundary conditions. 

Although this is an extremely simple model, the so- 
lution space is nevertheless surprisingly rich. We find gen- 
eralized paraboloidal solutions (§ 3.3), cylindrical solutions 
(§ 3.5), conical solutions (§ 3.6), and even converging solu- 
tions (§ 3.4). Moreover, some of the solutions have field lines 
crossing the Alfven surface, while others live entirely inside 
the Alfven surface. Numerical examples of the different kinds 
of solutions are shown in Figures 1-6, and the regions of pa- 
rameter space over which the different kinds are found are 
indicated in Figures 7 and 8. In §§ 3.3.1, 3.3.2 we discuss 
scaling relations for paraboloidal solutions, and in Figures 
9 and 10 we show that these analytical scalings agree well 
with numerical solutions. 

Given the bewildering variety of solutions, we are faced 
with the question: which is the "correct" solution for a given 
choice of v and Ml In an attempt to make some headway on 
this question, we may wish to impose the condition that the 
minimum fluid velocity i> m in, as defined in § 2.4, should not 
exceed c anywhere within the flow. However, this condition 
only eliminates half the parameter space of solutions (see 
Figs. 7, 8). Requiring in addition that the solution should 
be analytic on the axis, we find this corresponds to the de- 
sirable condition v m i n — > c asymptotically far from the disk 
(§ 3.3.3). This leads us to the stronger condition n m in = c at 
infinity instead of just ^ c, and it gives us a unique solution 
for all v ^ 1. However, for v < 1, it still leaves us with a 
continuous family of physically valid solutions which spans, 
for each choice of v and M, a range of values of s. To break 
this residual degeneracy, we are compelled to invoke yet an- 
other condition, viz., that the correct solution is the one 
with the minimum toroidal field, i.e., the minimum torque 
(Michel 1969). This latter condition subsumes the u m in — > c 
condition, and gives a unique solution for all values of v. 
It is, however, less well motivated than the analyticity con- 
dition on the axis or the requirement that the fast surface 
be located at infinity. We are thus led to use axisymmetric 
numerical simulations to verify which solution is picked out 
by the time-dependent system. 

§ 5 describes these numerical simulations and summa- 
rizes the results (see Figs. 11-15). For each of the three rep- 
resentative combinations of \y, M} that we analyze, we find 
that the numerical models converge to a unique final steady 
state, which is independent of the initial starting state. In 
all cases, the final state has a value of s equal to the smallest 
value of this parameter for which our analytical model gives 
a physically consistent solution (with v m i n ^ c over the do- 
main of interest) , i.e. , the solution with the minimum torque. 
This is also the solution with the largest acceleration. 

One interesting question is prompted by the numeri- 
cal simulations. In the mathematical analysis described in 
this paper and also in the physical discussion of the previ- 
ous paragraphs, we assigned great importance to the point 
at infinity, or at the edge of the grid, and assumed that 
one or more important boundary condition is set there. 
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This is a natural consequence of working with the steady 
state problem which involves analyzing the elliptic Grad- 
Shafranov equation. The numerical simulations, however, 
work with the hyperbolic time-dependent force-free equa- 
tions. Surprisingly, in the numerical simulations, the flow 
chooses the more-or-less correct configuration for the mag- 
netic field near the equatorial plane long before any signals 
(moving at the speed of light) have reached the outer bound- 
ary of the grid, or even reached the Alfven surface. Thus, the 
time- dependent problem appears not to require the point at 
infinity to set up the local flow. Exactly how does it do this? 
One answer is that the condition v m i u — > c at infinity is 
equivalent to analyticity on the axis, i.e., the absence of a 
singular current on the axis (§ 3.3.3). The numerical code 
forbids a singular current on the axis (at all z) and would 
therefore naturally pick out such solutions. For v < 1, we 
require an additional condition that the solution has the 
minimum torque (or maximum acceleration). This again is 
a condition that could, in principle, be applied at the equa- 
torial plane and perhaps this is how the code is able to pick 
the correct solution promptly. 

Ultimately, we suspect that only a perturbation analysis 
of the self-similar steady solutions we have obtained in this 
paper will provide real understanding of which solution is 
correct for each situation. It could be that the minimum- 
torque/maximum-acceleration solution is the only one that 
is stable to axisymmetric perturbations. 

The quasi-analytic work we presented is closely related 
to that of Contopoulos (f995a), who investigated self-similar 
force-free flows in order to better understand MHD solutions 
obtained by Li et al. (1992) and Contopoulos (1994). Our 
work differs from Contopoulos (1995a) in at least two ways: 
(i) We have identified a larger set of solutions ; (ii) Con- 
topoulos (1995a) suggests that the value of —sM (Ho in his 
notation) has to be fine-tuned in order to obtain solutions 
that reach asymptotic infinity with high Lorentz factors. For 
a physical, time-dependent solution, he suggests that this 
fine-tuning corresponds to the requirement that waves prop- 
agate from asymptotic infinity and reach the origin, and this 
would take an infinite amount of time. He suggests, there- 
fore, that such solutions are inapplicable as physical models. 

On the other hand, following Blandford (1976) for v — 
1, we have been able to identify a regularity condition on 
the polar axis that leads to a unique solution for each v ^ 1, 
and we use the minimum-torque (minimum-energy) condi- 
tion to choose a unique solution for each v < 1. As discussed 
above, we confirmed these conditions as natural by perform- 
ing time-dependent numerical simulations that generate so- 
lutions that are, by construction, regular on the polar axis 
and necessarily stable. These simulations indicate that the 
minimum-torque condition is the natural condition and that 
there is no fine-tuning required to obtain solutions that reach 
asymptotic infinity with w m i n — > c. However, for M close to 
unity the solutions with v < 1 are indeed sensitive to the 
value of s (Fig. 8). As applied to astrophysical winds, this 
may indicate that for disks with v < 1 the force-free approx- 
imation is difficult to maintain at arbitrarily large distances. 

Another important issue is how the self-similar force- 
free solutions described in this paper are related to self- 
similar MHD solutions such as those discussed by Li et al. 
(1992) and Vlahakis & Konigl (2003). We anticipate that 
MHD outflows that start out with high values of the magne- 



tization parameter a will closely follow one of our force-free 
solutions out to a certain distance before plasma inertia in- 
duces deviations. Exactly which of our force-free solutions is 
picked by the MHD flow? Following the work of Goldreich & 
Julian (1970), one suspects that MHD would automatically 
single out the minimum torque force- free solution, but this 
remains to be seen. It is also interesting to ask exactly where 
the MHD solution would start deviating from the force-free 
solution. 

Returning to the results obtained in this paper, we note 
that, to a good approximation, the thick solid line in Fig- 
ure 7 gives the minimum torque solution for any v and M. 
The eigenvalue s is given approximately by s ~ 2/u, and 
the magnetic field components and the Lorentz factor have 
the approximate asymptotic scalings given in §§ 3.3.1, 3.3.2. 
Some of these results have been derived by other authors 
(e.g., Beskin & Nokhrina 2006), but some appear to be new. 
The unification of v ^ 1 models and v < 1 models is another 
relevant contribution of the present paper. Our recent work 
(McKinney & Narayan 1995a,b) has shown that GRMHD 
models have currents and fields described by v « 3/4 for a 
wide range of conditions. It is therefore important to under- 
stand the properties of self-similar models with this partic- 
ular value of the index, originally introduced by Blandford 
& Payne (1982). 

An interesting result of our analysis is that rotation has 
almost no effect on the poloidal structure of field lines. This 
was known to be the case for the split monopole solution of 
Michel (1973a) and the v = 1 solution of Blandford (1976), 
but we now find that it is true for the entire family of self- 
similar solutions considered in this paper. It suggests that 
the collimation of astrophysical jets is not the result of the 
toroidal field associated with rotation. In our solutions, col- 
limation seems to be produced by the poloidal field itself. In 
effect, each field line is collimated by the pressure associated 
with field lines further out. However, this result is for the 
specific self-similar model we have considered, which has a 
flat rotation curve. It remains to be seen if the results carry 
over to a disk with a Keplerian rotation profile. 

In contrast to the case of collimation, our models show 
that rotation and toroidal field are critical for accelerating 
the force-free wind. Serious acceleration begins only when 
the toroidal component of the magnetic field dominates over 
the poloidal component, which happens only after a field 
line crosses the Alfven surface. The larger the rotation of 
the disk, the closer the Alfven surface is to the foot-point of 
a field line, and the sooner strong acceleration is initiated. 
This is illustrated in Figures 9 and 10. 

Finally, we note that the force-free winds we have con- 
sidered in this paper are highly idealized, and their relevance 
to real disk winds is unclear since MHD turbul ence within 
the disk may tangle up such large -scale fields ( Mc KinnevI 
120051: iMcKinnev fc NaravarJl2006bl) . Our hope is that some 
of the analytical results and qualitative insights obtained 
here may carry over to more realistic MHD models of winds. 
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APPENDIX A: ANALYSIS OF PARABOLOID AL 
SOLUTIONS 

In this Appendix we focus on paraboloidal solutions and 
derive a few asymptotic results in the limit of large it = 
z/R. We begin by separating out equation 12611 into terms 
of different orders: 

uTT" - u 2 T' 2 /is + (3-2v)uTT' + v(v - - s 2 )T 2 
- M~ 2 T {u+2)/v [u 2 T" + (3 - 2v)uT' + v(y - 2)T] 



+ TT"-T' 2 /u 



(Al) 



By "paraboloidal" solutions, we mean those in which T(u) 
goes as u~ M with a positive value of /i as u — » oo. We will 
concentrate on the region well outside the Alfven surface, 
where T(u) <C M v . It is obvious that the dominant terms 
in equation 1A11 in this region are those listed in the first 
line, and so we begin with an analysis of these terms. The 
single term in the fourth line is always much smaller than 
the others and we will neglect it. The terms in lines 2 and 
3 are of intermediate importance and we consider them as 
needed. 

The first line of equation ijAl^ is homogeneous and is 
identical to equation 153H in the main text. It has two in- 
dependent power-law solutions. We require a decaying solu- 
tion, and such a solution exists only when s > 1 (as in the 
main text, we restrict our attention to positive values of M 
and s). The solution is given by 



T{u) 



fj, = v(s — 1). 



(A2) 



In the following two subsections, we consider the effect of 
small perturbations on this solution and calculate the next 
order term in a power series expansion. 



Al Perturbation Analysis of the Homogeneous 
Equation 

Focusing still on just the homogeneous equation represented 
by the first line of l)Aljl . consider small perturbations around 
the power- law solution IA2I : 

T(u) w + e(u), (j, = v(s-l), u>u c . (A3) 

We will assume that e(u) <C Substituting (IA3I > in 1A1II 
and retaining only the leading terms, we obtain 

= u V'+f^+3-2i;W 



+ [v(v + l)-(3-2v)v + 2v(v-l)(l-s 2 )}e. (A4) 

This equation has power-law solutions of the form e ~ u~ , 
with two values of £: 



6 = n, 6 = H - - 1 + 2(1 - u) 



(A5) 



The first solution is trivial — it merely reproduces the 
leading order term in equation HA3t and carries no new 



information. The second solution is, however, interesting. 
It gives an asymptotic behaviour for T(u) of the form 



l (l + Cl 



), where C is a constant and 



£2 — M = 2(1 — u) I — + 1 



(A6) 



We now consider the stability of the solution; we mean 
stability in a mathematical sense, not dynamical. The per- 
turbation e(u) will remain smaller than the primary term 
u _M with increasing u only if £2 — A* > 0, i.e., only if v < 1. 
Thus, so long as v < 1, the paraboloidal solution is "sta- 
ble" for any positive value of /i, i.e., for any s > 1. This 
is confirmed by Figure 5, where we see that all the solu- 
tions with s > 1 are well-behaved and have power-law, i.e., 
paraboloidal, behaviour as u — > 00. 

When v > 1, however, we see that £2 — M is always 
negative, and so perturbations grow with increasing u and 
ultimately dominate over the leading order term u -M . This 
is a sign that the solution is "unstable," as confirmed by 
Figure 3. Apart from the cylindrical solutions, we see that 
nearly all the other solutions are either conical, i.e., T(u) 
goes to at a finite u = Uf, or singular, i.e., T(u) attempts, 
and fails, to recross the Alfven surface. Obviously, none of 
these solutions is physically satisfactory. There is, however, 
one (and only one) paraboloidal solution which manages to 
avoid the instability; it is shown as the thick solid line in 
Figure 3. To understand how this special solution manages 
to avoid growing perturbations, we note that perturbations 
to the primary solution « _M are generated by the terms in 
lines 2-4 of equation 1A1I . So let us briefly take a look at 
these terms. 

Consider the terms in the second line of equation 1A1I . 
They add up to zero when T(u) has the power-law form 
y -(2-i/) -pjjyg^ f or this particular power-law, these terms do 
not introduce any perturbations at all to the solution. In 
general, the primary power-law solution it -M is not of this 
special form. However, there is one particular case when 
the same power-law solution satisfies both lines 1 and 2 of 
equation HAlll . This is when 



2 

V 



/'• 



(A7) 



So for this one value of s, line 2 does not introduce any 
perturbations to the primary tt -M dependence of T(u). We 
would therefore expect the solution for this particular value 
of s to be stable. It is gratifying to see that the special 
solution shown by the thick solid line in Figure 3 does in- 
deed correspond to s being very nearly equal to this spe- 
cial value, and the solution does behave asymptotically as 
u _ ( 2_ "). Of course, we only considered the first two lines 
of lAH . For the above special solution, line 3 is of a lower 
order than line 2, but it is not completely negligible. There- 
fore, the paraboloidal solution does not correspond exactly 
to the value of s given in (IA7L but it is very close. 

In summary, for v < 1, a paraboloidal solution is avail- 
able for all s > 1, while for v > 1 it is available for only one 
specific value of s which is nearly equal to 2jv. 

A2 Power Series Expansion and Lorentz Factor 

In this subsection, we ignore line 2 of equation HAlll . When 
v > 1, we showed above that the contribution of this line 
vanishes for the only paraboloidal solution available. When 
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v < 1, this line is sub-dominant relative to line 3 for most 
paraboloidal solutions of interest to us, especially the ones 
on and above the thick solid line in Figure 7. 

Consider, therefore, the differential equation consisting 
of lines 1 and 3 of equation 1A11 , As before, consider a so- 
lution of the form 1A3II . where we assume that the second 
term is much smaller than the first. When this solution is 
substituted in line 1, the zeroth order terms cancel and the 
first order terms take the form given in the left hand side of 
equation 1A4L On the other hand, the zeroth-order terms 
give a non-vanishing contribution when substituted in line 
3. Thus, collecting terms, we obtain the following differential 
equation for e: 

2 a , / 2/x \ , 

it e + h 3 — 2v \ue 



+ [nin + 1) - (3 - 2v) l i + 2v{v - 1)(1 - s 2 )} 



Solving this, we obtain 
e(it) = -Cm- (m+2) , 



-fi-2 



iu, = v(s-l). (A8) 



(A9) 



where C is a coefficient which is expressed in terms of v, fi 
and s. In the rest of this subsection we will ignore coefficients 
of order unity, such as C, since our primary interest is in the 
scalings of various quantities. Thus, we write the solution for 
T(u) as 



T(u) 



(A10) 



Clearly, this solution gives the first two non-zero terms of 
a power series expansion in 1/it of the paraboloidal solu- 
tion (note that the term proportional to 1/tt vanishes), i.e., 
the first two nonvanishing terms of an expansion around the 
point at infinity. For many purposes, the first term is good 
enough, but the second term is essential for estimating the 
Lorentz factor of the force- free outflow, as we now discuss. 
Incidentally, to get an idea of the error we make by neglect- 
ing coefficients of order unity, note that the power series 
expansions of the two analytic solutions given in §§ 3.1 and 
3.2 are given by 



T(u) 



t(„) « v = ' " = % (A12) 

The factors 1/2, 1/4 and 3/4 are omitted when we write the 
series expansion approximately as in ijAlO^ . 

The stream function P(R, z) corresponding to the above 
solution is 



P(R,z) = R u T(u) 



R' 



(A13) 



Since P is constant on field lines, we can solve for the shape 
of a field line with footpoint at R{ p : 



R 

Rt P 



z 

R~Tp 



(s-l)/s 



1 + 



R fp 

z 

Rt P 

z 



2v/(p, + v) 



2/s 



(A14) 



where, in the final expression, we have used equation 1A21 
to substitute for /j,. From equations 127L 128L 129L we can 
estimate the three components of the magnetic field: 



B R (R,z) 
B^R,z) 
B z (R,z) 



u(v 2 —2v + fiv — n)/i' 

—M— r ^- 7 

zv \ z 2 



R 



(A15) 
(A16) 
(A17) 



We have written the first two terms for and B z , but only 
the leading term for Br since this component of the field is 
very weak at large u. 

Consider now the variation of the field components as 
a function of z along a given field line. To calculate this, we 
substitute for R from equation l|A140 : 



Mi) 



B s 



\R tp 



Mi) 



z 
M 

Rip 



(2s-l)/s 



-Rfp 

z 



(s-l)/s 



2/s 



2(s-l)/s 



-Rfp 



2/6 



(A18) 



(A19) 



(A20) 



From these, we obtain the following expressions for the 
poloidal and total field strength along a field line: 



Bl {~k 



Rip 

z 



4(s-l)/s 



Rip 

z 



2/s 



(A21) 



B' 



(.Rip) 



M 



R 



2(s-l)/s 



(A22) 



R, 



2/s 



M 2 \ z 



2(s-l)/s 



Substituting these expressions and (I A 1411 in equation 13611 . 
we obtain the following estimate for the minimum Lorentz 
factor along a field line, 



In 



2/s 



+ 



M 2 



Rap 

z 



2(s-l)/s 



-1/2 



(A23) 



Since there are two terms, we need to keep track of the 
behaviour of each. 

First, consider the case v > 1. In this case, we saw 
earlier that 1/s = u/2 > 1/2. Therefore, for large z/R{ p , the 
second term in equation l|A23fl always dominates, and so we 
have 



M 



(s-l)/s 



(2-v)/2 



v > 1, (A24) 



-Rfp / V B-fp / 

where we have substituted for s using equation ljA7fl . For 
v — 5/4, this gives 7 m i n oc z 3//8 , which agrees very well 
with the scaling we find for the numerical solution shown in 
Figure 4 (thick solid line). 

Consider next v < 1. In this case, the second term in 
equation 1A23> dominates for a range of z/R{ p up to a cer- 
tain limit and the first term dominates beyond that, i.e., 



M 



-Rii 



(s-l)/s 



U<1, 



z < 2 crit ,(A25) 
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7min « \R~) ' I/<1 ' Z>2 cr it, (A26) 

where the transition value of z is given by 

2crit =R^M- s,{s - 2) . (A27) 

The two regimes of 7 m in are seen clearly in the bottom right 
panel of Figure 9. 

Before concluding, we remind the reader once again that 
all the relations given in this subsection have the correct 
scalings, but we have ignored numerical coefficients. For in- 
stance, the coefficient M in the expressions for B$ and 7mm 
is perhaps more appropriately written as sM. However, s 
is a quantity of order unity and so, in the general spirit of 
neglecting coefficients, we have ignored this refinement. The 
other point is that we ignored line 2 of the differential equa- 
tion lAH in the discussion here, whereas the terms in this 
line are negligible only for a limited range of s. A more com- 
plete analysis would retain these terms, but this is beyond 
the scope of the paper. 
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